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Abstract 


Geometric and Material Nonlinearities which arise due to large deformations are mod- 
elled in the present work using finite element method. Modelling of large deformations 
is based on Total Lagrangian formulation. Materials having elastic-perfectly plas- 
tic and elastic-linearly work hardening characteristics are considered in the present 
study. Fracture parameters, Energy Release Rate and Crack Tip Opening Displace- 
ment (CTOD) are evaluated for a simply supported beam with central crack and a 
plate with center and edge crack of finite dimensions. 
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Chapter 1 


INTRODUCTION 


1.1 Importance of the Problem 


Fracture mechanics supposes the pre-existence of a significant crack like defect that 
will lead to failure. Elastic or eleistic plastic continuum mechanics theory is used 
to develop parameters which characterise the crack tip stress field. If a material 
property is measured in terms of suitable parameters under known circumstances, 
fracture mechanics then allows this property to be used to calculate the behavior of 
other geometrical configurations for similar physical circumstances. 

When a body having cracks is subjected to loading there will be large deformation 
and also a plastic zone in front of the crack tip. Much work has been done in recent 
years to estimate the fracture parameter of interest in cracked bodies. However, 
most of the earlier analyses were based on linear elastic behavior and assume the 
deformations to be small. Some work was also carried out to study the effect of large 
deformations on fracture parameters. However, most of this was limited to Linear 
Elastic Fracture Mechanics (LEFM). 

Unfortunately, in many engineering problems it is indispensable to consider the 
effects of geometric and material nonlinearities. In metal-forming processes, work- 
pieces may undergo quite large displacements, large rotations, and large deformation. 
In addition when the stress or strain within the material reaches a certain level, 
the deformation of the material will no longer be linearly elastic, and non-linearly 
elastic or plastic deformation wiU occur, with plastic deformation, the stress-strain 
relation is no longer unique and will depend on deformation history. In fact, even 
when the global deformation of a structure can be regarded as elastic, local plastic 
deformation can often be found in the structure. This is especially true when con- 
tacts or cracks are present, since contacts can often produce high local pressures on 
contacting boundaries and there will be a plastic zone ahead of the crack tip. More 
importantly, many contact-impact problems of great interest involve plastic deforma- 
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tion. In metal-forming processes, for example, plastic deformation of the work piece is 
involved without exception. Another important case in which plastic deformation is 
involved is structural crashworthiness analysis in which structures are often subjected 
to heavy loadings. 

Thus, it can be concluded that it is very important to model large elastic-plastic 
deformations in engineering problems. The present work is aimed at modelling large 
elastic-plastic deformations and estimating the effect of geometric nonlinearity on 
fracture parameters for two dimensional bodies under plane stress/plane strain con- 
ditions. 


1.2 Literature Survey 


Lagrangian formulation is by far the most widely used formulation in finite element 
applications to nonlinear structural problems [1-6]. Many of the early incremental 
formulations are incomplete in that the variations in loading and rotations of the 
element are not properly represented [5]. The most general approach to incremental 
lagrangian formulation is to consider a fully nonlinear kinematic relations within a 
linear increment [5,6]. 

When fully nonlinear kinematic relations are used within a linear increment, the 
number of stiffness contributions to the nonlinear element stiffness matrix are found 
to be three matrices [5,6]. These are the usual small displacement or the incremental 
stifi'ness matrix, the initial stress or the geometric or the tangent stiffness matrix and, 
finally the initial displacement or the initial rotation stiffness matrix. 

Geometric nonlinearities were first included by means of incremental geometric 
stiffness [7-9]. Later two different approaches were pursued in incremental nonlinear 
finite element analysis. In the first approach referred to as the moving coordinate 
formulation, the static and kinematic variables were referred to the deformed config- 
uration in each load step [11,12]. [13] gives some results based on this approach. 

In the second approach, termed as the lagrangian formulation, the static and 
kinematic variables were referred to the initial configuration of the body. [5,6] used 
this approach. 

Bathe et al [15] presented a more general finite element formulation based on 
continuum mechanics principles accounting for all nonlinearities and showed that 
the finite element solution based on the updated and total lagrangian configurations 
remained the same, provided proper constitutive relations are employed. It is the 
numerical effectiveness which governs the choice between the total and updated la- 
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grangian configurations. 

In the study of fracture mechanics interest is often focussed on the singularity point 
where such quantities as stress become mathematically infinite. Crack tip elements 
were designed to simulate the singularity at the crack tip. [16-19] studied the use of 
crack tip elements in fracture problems. It was reported that the crack tip singularity 
characteristic of LEFM can be obtained in a 2D 8-noded isoparametric elements by 
shifting mid side nodes to the quarter point. The use of finite elements in fracture 
mechanics increased. But, many of the results were based on LEFM and do not 
consider the nonlinearities present in the structure. 

[20] determined the combined effect of shear and large deformation on stress in- 
tensity factors. The influence of geometric nonlinearity on fracture parameter was 
found to be considerable. However, the study was still limited to LEFM. 

In Elastic Plastic Fracture Mechanics(EPFM) J-Integral and CTOD are the gen- 
erally used fracture parameters. [21] finds the J-criteria to be most attractive from 
computational viewpoint. [22] discusses the approaches to establish fracture criteria 
which would allow the engineer to decide whether a structure is safe against steady or 
unsteady crack growth under given conditions. [23] shows that calculation results ob- 
tained with the method of finite elements, employing an increrhental plasticity model 
can serve to evaluate the independence of J. Crack tip opening displacement was also 
favored by some researchers. [24] studies the effect of thickness on crack opening dis- 
placement. [25] carries out fractographic studies of crack tips in steel and aluminium 
alloys. 

[26] describes how triangular quarter point elements can be used as elastic as well 
as perfectly plastic crack tip elements. The elastic singularity (l/>/r) is obtained by 
having one node at the crack tip. The plastic singularity (l/r) is obtained by having 
multiple independent nodes at the crack tip. Geometric Nonlinearity is modelled 
using UL formulation in [35]. 


1.3 Outline of the Thesis 


Importance of the problem and literature survey are discussed in chapter 1. Chapter 
2 summarises the formulation of the general incremental equations of motion starting 
from continuum mechanics principle. The methodology for elastic-plastic analysis is 
explained. The total lagrangian formulation on which the present work is based is 
also discussed in detail in this chapter. Finite element implementation is discussed 
in chapter 3. Crack tip singularity modelling and evaluation of fracture parameters 
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is explained in chapter 4. Comparison of our results with available analytical results 
is given in chapter 5. Some more sample problems are also solved in this chapter. 
Chapter 6 gives conclusions and scope for further improvement in the present work. 



Chapter 2 


FORMULATION OF THE 
PROBLEM 


The objective of this chapter is to give a continuum mechanics based incremental 
equations of motion which form the basis for general nonlinear finite element analysis. 
The equations of motion are based on the principle of virtual work which is dealt with 
in Section 2.1.1. The stress and strain tensors which are found to be appropriate in 
the incremental formulation of the virtual work principle are introduced in section 
2.1.2. Section 2.1.3. deals with total lagrangian formulation which is widely used in 
large deformation analysis. The methodology for elastic-plastic analysis is given in 
section 2.1.4. 


2.1 Formulation of the Continuum Mechanics In- 
cremental Equations of Motion 


In the conventional linear analysis of a body, the deformation of the body is assumed 
to be small, which means, the equations of motion for any applied load can be set up 
in the initial configuration of the body. But, in the more general analysis involving 
finite deformations, the equilibrium of the body has to be established in the cur- 
rent configuration. Unfortunately, the current configuration is what is the unknown 
quantity. This calls for the incremental analysis to be performed in the case of large 
deformation problems. These sections in this chapter aim at deriving the incremental 
equations of motion based on the virtual work principle. 
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2.1.1 Principle of Virtual Work 

Consider the motion of a body in a cartesian coordinate system as shown in Fig. 
2.1. It is assumed that the body can experience large displacements, small strains, 
large rotations and a nonlinear constitutive response. The aim is to evaluate the 
equihbrium states of the body at discrete times 0, At, 2At, 3At, where At is an 
increment in time. 

The equihbrium of the body at time t + At based on the principle of virtual 
displacements is given by 

/+Aty = (2.1) 

In the above equation and the equations to follow, the left hand superscript de- 
notes the configuration at which the quantity is measured and the left hand subscript 
denotes the configuration which is taken as reference for measurement. In case the 
quantity considered occurs in the same configuration in which it is measured, the left 
subscript is not used. Thus, are the cartesian components of the cauchy stress 

tensor at t + At and t+At^ij ‘'■re the cartesian components of infinitesimal strain tensor 
due to the virtual displacements 6u and S means ’’variation in” i.e., 

- _ A ( dui duj N 1 / dSui dSuj 

2 ^ d*+^*Xi 

It should be noted that the cauchy stresses are ’’body internal forces per unit 
area” in the configuration at time t + At and the infinitesimal strain components are 
also referred to the equilibrium configuration at time t + At. L.H.S. of virtual work 
expression (2.1) is the internal virtual work performed when the body is subjected 
to virtual displacements at time t + At and the corresponding external work due to 
virtual displacements is given by 

6ui Sui ( 2 . 2 ) 

where and fi are the components of externally applied body and surface 

force vectors at time t + At with reference to the same configuration. It is to be noted 
that the summation convention of tensor notation is followed throughout. 

The difficulty in solving equation (2.1) is that the configuration at time t -f At is 
not known. This is an important difference compared with the linear analysis where 
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it is assumed that the displacements are so small that the configuration change is 
negligible. 

The continuous change in the configuration of the body with loading in a large 
deformation analysis is dealt with in an elegant manner by using appropriate stress 
and strain measures, as discussed in the sections to follow. 


2.1.2 Stress and Strain Tensors 

As was pointed earlier, in the nonlinear analysis, special attention must be given to 
the fact that the configuration of the body changes continuously. This is taken care by 
choosing proper stress and strain measures. The stress and strain measures chosen 
should express the internal virtual work in equation (2.1) in terms of an integral 
over a volume that is known and it should be easy to incrementally decompose the 
stresses and strains in an effective manner. Though there are various stress and strain 
measures available, in principle, the stress and strain tensors that will be most suited 
for incremental formulation of virtual work principle are Second Piola-Kirchhoff stress 
tensor and Green-Lagrange strain tensor [15]. 

The Second Piola-Kirchhoff stress at time t referred to the configuration at time 
0 is given by 


t 

0 



V 




0 

t 


Xj,n 


(2.3) 


where is the derivative of the coordinate in configuration at time 0 with 

respect to coordinate and 7^ represents the ratio of mass densities at time 0 and at 
time t. The determination of the ratio of mass densities is explained in the appendix. 
From equation (2.3) it can be seen that, by a purely kinematic transformation, the 
second Piola-Kirchhoff stress tensor at any configuration can be calculated from the 
corresponding cauchy stress tensor and vice versa. 


The strain tensor used with the second Piola-Kirchhoff stress tensor in the large 
deformation analysis is the Green-Lagrange strain tensor and it is defined at any time 
t as 


1 

in which lujj = is the derivative of the displacement component in configuration 
at time t with respect to coordinate 
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The stress and strain measures described are symmetric and their components are 
unaffected by rigid body rotations of the material [27]. A most important aspect is 
that the stress and strain tensors described above are energetically conjugate to each 
other [27], which is a necessary condition for using them in incremental formulations. 


2.1.3 Total Lagrangian Formulation 

Consider the deformation of the body starting from the configuration at time 0 to 
the configuration at time t + At passing through discrete time steps At, 2At , ... as 
shown in the Fig. (2.1). In developing the solution for the virtual work expression 
(2.1) established in the configuration at time t + At, it is assumed that the solutions 
for the static and kinematic variables for all time steps starting from time 0 to time 
t, inclusive, have been obtained. The aim is to evaluate the solution variables at time 
t + At. After achieving this, the same strategy can then be applied repetetively till 
the complete solution path is solved. Hence, in the analysis the path of the body is 
followed from the original configuration (i.e. at f = 0) to the final configuration of 
the body. This method of approach is referred to as the lagrangian formulation [15] 
and widely used in the nonlinear analysis. 

It was discussed in section 2.1.1 that the solution to virtual work expression (2.1) 
can be obtained if all the variables in the expression can be referred to some known 
configuration. For this purpose, in principle, any one of the already calculated equi- 
librium configurations could be used. However, the choice lies between two different 
formulations, namely, the total lagrangian formulation (T.L) and the updated La- 
grangian formulation (U.L). In the total lagrangian formulation the solution variables 
in the configuration configuration t + At can be obtained by referring all the variables 
to the undeformed equihbrium configuration at time 0. In the updated lagrangian 
formulation, the solution at time t At can be obtained from the configuration at 
time t i.e., the variables are all referred to the configuration at the end of the previous 
time step, which represents the load step in the case of static analysis. 

The virtual work expression (2.1) can be transformed in terms of second Piola- 
Kirchhoff stress and Green- Lagrange strain as 

= (2.5) 

where are the cartesian components of the second Piola-Kirchhoff stress 

tensor and are the cartesian components of the Green-Lagrange strain tensor 

at time t + At referred to the configuration at time 0. In the analysis the load is 
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assumed to be deformation independent and hence is calculated using 

Sui ^+^^clA + [ Sui (2.6) 

Jt+At^ Jt+Aty ^ ' 

Approximate solutions to equations (2.5) and (2.6) can be obtained by linearizing 
the relations which is explained below. 

1. Equation of Motion 



Joy = 

(2.7) 

where 

<+A*C 0 t+At 0 

0 — t+At^ t-hAt'^^yrn f mn t^At'^ 3 ,n 

r 

(2.8) 



(2.9) 

Incremental Decompositions 
a) Stresses 


b) Strains 

t-f At c i C _L Q 

0 ‘ 

(2.10) 


. _L 1 

0 1 

= o^tj “1" o^z'i J 

(2.11) 


O^ij — “b O’^i.t "b 0^k,i O^ArJ "b 1 /2 22) 

07»i ~ 2^'^k,i O'^kj'i J 

This can be derived using taylor’s series 
3. Equation of motion with Incremental decompositions 

= <Soeii 

oSiy Body Uv + ISii SoViJ °dv = ‘+^‘R - ISyy So^y °dV (2.13) 
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4. Linearisation of Equations of Motion 
using the approximations 

oSij — O^^ijrs O^rsj ^O^ij ~ 

we obtain the approximate equation of motion: 

( 2 . 14 ) 


jQy O^ijrs O^rs ^O^ij ^O^ij 

ISij °dV 


where 


Q^ijrs — ^ ^ij ^Ts + ^js + ^is ^jr ) 

where A and (x are Lame’s constants and % is the kronecker delta. 

E V E 


A = 


(1 + 1^) (1 — 2i/)’ ^ 2(1 + 1/) 


The relation (2.14) can be employed to calculate an increment in the dispacements, 
which then is used to evaluate approximations to the displacements, strains and 
stresses corresponding to t + At. The displacement approximations corresponding to 
t + At are simply obtained by adding the calculated increments to the displacements 
of time t, and the strain approximations are evaluated from the displacements using 
the available kinematic relations i.e 

Q^ij = 2^o'^hj "k o%>' "k 

The stress can then be calculated using the constitutive relations for small strain 

i Q, , t i , 

Q^ij — Q^ijrs o^rs 

Assuming that the approximate displacements, strains and thus stresses have been 
obtained we can now check into how much difference there is between the internal 
virtual work when evaluated with the calculated static and kinematic variables for 
time t+At and the external work. Denoting the approximate values with a superscript 
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(1) in anticipation that an iteration will in general be necessary, the error due to 
linearisation is 


Error °dV (2.15) 

The above considerations show that the R.H.S. in (2.14) represent an ” out-of- 
balance virtual work ” prior to the calculation of the increments in the displacements, 
where as the R.H.S. of (2.15) represent the ” out-of-balance virtual work ” after the 
solution, as the result of the linearisations performed. In order to further reduce the 
” out-of-balance virtual work ” we need to perform an iteration in which the above 
solution step is repeated until the difference between the external virtual work and 
the internal virtual work is negligible within a certain convergence measure. Using 
the T.L. formulation the equation solved repetitively, for k = 1,2,3, .. . is 

f oCijrs Aoe("; doe,-,- °dV + f ^Sij = 

(2.16) 

where the case ^ = 1 corresponds to the relation (2.14) and the displacements are 
updated as follows 

«+At„(0) ^ 

Equations above correspond to modified newton-iteration scheme. 


(2.17) 


2.1.4 Methodology for Elasto-Plastic Aii 2 dysis 

Before the onset of plastic yielding the relationship between the stress and strain is 
given by the standard linear elastic expression 

(Tij = Cijki tki (2.18) 

where (Tij and tki are the Second Piola-Kirchhoff stress and Green Lagrange strain 
components respectively and Cijki is the tensor of elastic constants which for an 
isotropic material has the form 
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^ijkl — ^ ^ij ^kl ”1“ ^ik ^jl “t" P ^il ^jkj 
where A and fj, are the lame constants and Sij is the kronecker delta. 

Plastic behavior is characterised by an irreversible straining which is not time 
dependent and which can only be sustained once a certain level of stress has been 
reached. 

To model elasto-plastic material deformation three requirements have to be met: 

• An explicit relationship between stress and strain under elastic conditions must 
be formulated. 

• A yield criterion indicating the stress and strain level at which plastic flow 
commenses must be presented. 

• A relationship between stress and strain must be developed for post yield behav- 
ior i.e when the deformation is made up of both elastic and plastic components. 

The relationship between stress and strain under elastic conditions is given by (2.18). 

Yield Criterion 

The yield criterion determines the stress level at which plastic deformation begins 
and can be written in the general form 

f(a,i) = K{k) (2,19) 

where / is some function and AT is a material parameter (which may be a function of 
the hardening parameter k) to be determined experimentally. On physical grounds, 
any yield criterion should be independent of the orientation of the coordinate system 
employed and therefore it should be a function of the three stress invariants only. 

Jl = (Tii 1 

J 2 = I (2.20) 

Experimental observations, notably by [28], indicate that plastic deformation of 
metals is essentially independent of hydrostatic pressure. Consequently the yield 
function can only be of the form 
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where J 3 ^'I'e the second and third invariants of the deviatoric stresses. 

There are many yield criteria mentioned in the literature, but Tresca and Von 
Mises criteria are very popular. 

Tresca Criteria 

According to this criteria yielding begins when the maximum shear stress reaches 
a certain value. If the principal stresses are (7i,<T2,<r3 where ctj > <t2 > aZ then 
yielding begins when 


<y\ — (7z = Y{k) (2.21) 

where V is a material parameter to be experimentally determined and which may be 
a function of the hardening parameter k. 

Von Mises Criteria 

Von Mises suggested that yielding occurs when J'^ reaches a critical value, 

( 4 )''“ = K{k) ( 2 . 22 ) 

in which K is a material parameter to be determined. can be explicitly written as 

= \ ^'ij = ^[(^1 - <^ 2 ? + (<72 - < 73 )^ + (<73 - <7i)^] (2.23) 

Yield criterion may be further written as 

W = ^/Z = VZ K (2.24) 

where a = ^ {crb is termed the effective stress. 

Elasto-Plastic stress strain relation 

After initial yielding the material behavior will be partly elastic and partly plastic. 
During any increment of stress, the changes of strain are assumed to be divisible into 
elastic and plastic components, so that 


dtij = + (deij)p 


(2.25) 
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where 


, , . 5cr' ■ (1 - 2z/) . , , , 

= ^ + ■ E 

In order to derive the relationship between the plastic strain component and the 
stress increment a further assumption on the material behavior must be made. It will 
be assumed that the plastic strain increment is proportional to the stress gradient of 
a quantity termed the plastic potential Q, so that 

(&«), = " ^ ( 2 - 27 ) 

where d\ is a proportionality constant termed the plastic multiplier [29]. The above 
equation is termed the flow rule since it governs the plastic flow after yielding. 

The identity f = Q gives us the normality condition[30] 

= (2.28) 

and is termed the normality condition since is a vector directed normal to the yield 
surface at the stress point under consideration. Experimental observations indicate 
that the normality condition is an acceptable assumption for metals. 

Substituting and (deij)^ in (2.25) gives 

= ^ + + (2.29) 

For matrix formulation the yield function f{ff) = K{k) can be rewritten as 

/(o) - K(k) = 0 


or 

Fia, k) = f{a) - K{k) (2.30) 


By differentiating we have 




dF 

dcr 


dcr + 


dF 

dk 


dk = 0, 


(2.31) 
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or 


where 


and 


dcr — A dX = 0 


T dF ,dF dF dF 


r 

« = ^ = [ 


dcr day’’ dr^y 


A = -— — dk 

^ iX dk 


Therefore, equation (2.29) reduces to 


(2.32) 


(2.33) 


(2.34) 


1 iJ T 

dt = [DY^da FdX— (2.35) 

where [D] is the matrix of elastic constants. Premultiplying both sides of (2.35) by 
do = aY D and eliminating da by use of (2.32) we obtain the plastic multiplier 


dX = T = r oY dn dt 

[AAaT D a] ^ 

Substituting this in (2.35) we get the complete elasto-plastic incremental stress strain 
relation to be 

da = Dep dt, 


where 


Dep = D — 


D a aY D 
A + D a 


It can be shown that [30] 


A = 



which is the local slope of uniaxial stress/plastic strain curve and can be determined 
experimentally from 


H' = 


Et 


1 _ Et. 

^ E 


where Et is the elastic-plastic tangent modulus of uniaxial stress-strain curve. 



Chapter 3 


FINITE ELEMENT 
IMPLEMENTATION 


In the previous chapter the incremental continuum mechanics equations that form 
the basis of general nonlinear displacement based finite element analysis were pre- 
sented. In this chapter we use these equations to develop the governing finite element 
equations. The numerical evaluation of elasto-plastic material matrix will also be ex- 
plained. At the end the numerical algorithms employed in solving resulting nonlinear 
finite element equations and the convergence criteria are dealt with. 


3.1 Isoparametric finite element discretisation 


By invoking the principle of virtual displacements for each of the nodal point displace- 
ments in turn, the governing finite element equations are obtained. The basic steps 
in the derivation of the governing finite element equations are the same as those used 
in linear analysis: the selection of the interpolation functions and the interpolation 
of the element coordinates and displacements with these functions in the governing 
continuum mechanics equations. 

In the incremental nonlinear analysis, since the new element coordinates are ob- 
tained by adding the element displacements to the original coordinates, it follows that 
the use of the same interpolations for the displacements and coordinates represent a 
consistent solution approach. 

Since the aim of the present work is to evaluate the effect of geometric nonlinear- 
ity on fracture parameters for two dimensional plane stress/plane strain conditions, 
attempt is made here to develop the finite element matrices for 8-noded isoparametric 
element with two degrees of freedom {u,v) at each node. 

Consider an 8-noded element in the configuration at time 0 and at time t in the 
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global coordinate system as shown in Fig. 3.1. In the isoparametric element solution, 
the coordinates and displacements are interpolated usi 

z = 1,2 

\ = h, 

k=l 

i = 1,2 

where is the coordinate of nodal point k corresponding to direction i at time t 
and is the displacement of nodal point k corresponding to direction i at time t. 
hk is the interpolation function corresponding to nodal point k and N is the number 
of element nodal points. 

Using equations (3.1) and (3.2) in equation (2.14) to evaluate the displacement 
derivatives required in the integrals, we get considering a single element 


(3.1) 


(3.2) 


(‘„Ki + IKnl) (Au<‘>) = •+^‘B - (3.3) 

where qKl is linear strain incremental stiffness matrix, 

qKnl is Nonlinear strain (geometric or initial stress or tangent) incremental stiff- 
ness matrix, 

t+Atp{x-i) -g Yg(.^Q;i. gf nodal point forces equivalent to the element stresses at time 
t -|-At and iteration (* — 1) and 

ig vector of externally applied nodal point loads at time t + At. 


Sifj, a = (for IBl oC '„Bl ’’dV) u 1 
IKnl a = (foy IBlt, iS IBn UV) a , 
IF = Soy iBl IS W J 


(3.4) 


qKi,, IKnl and IF are obtained from 


foy O^ijrs O^rs ^O^ij 

foy ‘oSii s„nii °dv 

I.y 'oSii Socij W J 


(3.5) 
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respectively. 3x3 gaussian integration is used in the present analysis to evaluate the 
linear and nonlinear stiffness matrices. 

The explicit form of the finite element matrices are given in the appendix. 


3.2 Numerical Computation of Elasto-plastic ma- 
terial matrix 


For numerical computations it is convinient to rewrite the yield criteria considered 
in chapter 2 in terms of alternative stress invariants. This formulation is due to [31] 
and its main advantage is that it permits the computer coding in a general form and 
necessiates only the specification of three constants for any individual criterion. 

We define a parameter 0 to be a convinient alternative to the third invariant Jg . 
The physical interpretation of $ is evident from Fig. 3.2. 

Based on these alternative stress invariants we can rewrite the Tresca and Von 
Mises yield criterion as 

2 cosO = crY{k) (3-6) 

and 

(3.7) 

respectively. 


In order to calculate the elasto-plastic material matrix we require to express the 
flow vector 'a' in a form suitable for numerical computation. It can be shown that 
the flow vector 

“ da 




dF dJi 


+ 


dF 


dJi da d{J^y/^ 


da d6 da 


can be written as [30] 


a — Cl Oi + C 2 02 T ^3 O 3 


(3.8) 


where 


= { 1 , 1 , 0 } 



Line of pure 
^hear {6 = 0) 


Figure 3.2; 


Two dimensional representations of the Tresca and Von Mises yield criteria 


CENTRAL library 

\ f. T.. KAWR 

111319 


iMLlfii. A 


• m •mm» 
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<^2 — 2J/1/2 

= {(y).( + y)}> 

I 

Cl = 0 C 2 = 2 cos0(l + ian^ tan36) C 3 = ^ (^^) tresca criterion. 

Cl = 0 C 2 = \/3 C 3 = 0 for Von Mises criterion. 

Thus the elasto-plastic material matrix can be determined from 


r>ep = O - 


D a D 
A + a'^ D a 


where A = H' the strain hardening parameter. 

Once the large deformation analysis is implemented and the elasto-plastic ma- 
terial matrix is calculated, assuming small strains large elastic plastic deformations 
can be modeled by simply substituting the Second Piola Kirchhoff stress tensor and 
Green-Lagrange strains in the material characterisation for the small displacement 
engineering stress and strain measures. 


3.3 Numerical Procedures for the solution of Non- 
linear finite element equations 


Having carried out the selection of appropriate finite element model, for the nonlinear 
analysis, a reliable and accurate response prediction of the model is essential to make 
use of the results with confidence. In the incremental solution procedure, it is of 
utmost importance to satisfy the governing finite element equations in each load 
increment to sufficient accuracy. Otherwise, the solution errors may accumulate and 
lead to solution instabilities. 

The basic equations to be solved in the nonlinear analysis corresponding to time 
t + At are 


‘+'^‘B(u*)-’+"F(u') = 0 


(3.9) 
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where is the vector of externally applied nodal point loads corresponding to 

time t + At and is the internal force vector in the same configuration. If the 

load increments are small the approximate equation to be solved reduces to 

But if a large increment is used the solution of equation (3.9) involves iteration. 

Two methods are widely used in solving nonlinear finite element analysis equa- 
tions: Newton- Raphson and the Modified Newton- Raphson method. The only dif- 
ference between these two methods is that in Newton-Raphson method the stiffness 
matrix is updated in every iteration, whereas Modified Newton-Raphson method 
employs lesser number of stiffness updations and hence may need more number of it- 
erations than the full Newton-Raphson method. The details about these two methods 
are explained in [27]. 


3.4 Convergence criteria 


In an incremental solution procedure, a realistic criteria to terminate the iteration 
process is important. After the solution variables are obtained after each iteration, a 
check has to be made to see whether the solution has converged within the specified 
tolerance. 

The tolerance can be specified on displacements, forces or stresses: at the end of 
each iteration a check is made to see whether the increment in displacement, or norm 
of out of balance load vector or increment in stresses is within specified tolerance 
respectively. 

In the present work the criterion for convergence is based on out of balance forces 
at the end of each iteration. This can be explained as 


'JzL 


X 100 < TOLER 


where N is the total number of nodal points in the problem and r denotes the iteration 
number. 


//s are the applied forces and are the residual forces. 
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This criterion states that convergence occurs if the norm of the residual forces 
becomes less than TOLER times the norm of the total applied force. The multipli- 
cation factor of 100 on the left hand side allows the specified tolerance factor to be 
considered as a percentage term. 

The tolerance should be specified in such a way that the results should be ob- 
tained with sufficient accuracy but at the same time the process should not become 
computationally expensive. 

3.5 Program Implementation 

This section deals with the computational steps involved in developing the software 
for modelling large elastic plastic deformations of two dimensional plane stress/plane 
strain problems. The algorithm is presented with reference to the Newton-Raphson 
and Modified Newton-Raphson technique. 

Computational Steps: Iteration number is represented by i and load increment 
by m. At start, m = 1 and z = 1. 

1. Reading the applied incremental force vector so that the total external force is 

given by + ™{AR}. 

2. Checking if the iterative technique is Modified Newton-Raphson method and 
m 1. If Yes, go to step 4, else the full Newton-Raphson method is employed 
by going to step 3. 

3. Calculating the overall stiffness matrix '^[K] = -|- '^\K]nl for the struc- 

ture using the equations developed in this chapter. The elements of linear and 
nonlinear strain displacement transformation matrices used in the determina- 
tion of the linear and nonlinear stiffness matrices are given in the appendix. The 
elastic or elasto-plastic material matrix has to be accordingly used depending 
on whether an individual element has yielded or not. 

4. Solving for the incremental displacements u in equation (3.3) where '^R is the 
force vector used in the R.H.S. 

5. Calculating the incremental displacements, strains and stress components. 

6. Calculating the total stresses at the end of iteration i of load increment m. 

7. Finding the equilibrium force vector using stress vector calculated in step 6. 
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8. Using the force criterion to check for convergence. 

9. If the solution has converged proceed with next load increment else find unbal- 
anced force vector and go to step 2. 

Stop when all load increments are finished. 



Chapter 4 

EVALUATION OF FRACTURE 
PARAMETERS 


Cracks are present to some degree in all structures. They may exist as basic defects 
in the constituent materials or they may be induced in construction or during service 
life. Therefore, to prevent failure by crack propogation it is necessary to evaluate 
fracture parameters. The fracture parameters considered in the present study are 
Energy release rate and CTOD. Analytical solutions exist for selected, relatively 
simple, cases, and it is in the determination of the stress and displacement fields 
for complex practical situations that numerical techniques such as finite elements 
and boundary integral methods play an important role [32]. Crack tip singularity 
modeling is discussed in section 4.1. Estimation of fracture parameters is discussed 
in section 4.2. 


4.1 Crack Tip Singulatrity Modeling 


LEFM studies have shown that the stresses nearer to the cracktip vary inversely with 
the square root of the radial distance (r) from the crack tip. This implies that the 
stresses at the crack tip posses 1 / yjr singularity. For a elastic perfectly plastic material 
once we come to the plastic region stresses at the crack tip posses 1/r singularity. For 
a material with strain hardening the singularity is of the type 1/r” where .5 < n < 1. 
So elements near the crack tip should be able to model 1 / \Jr singularity while within 
the elastic limit and 1/r singularity after crossing the elastic limit. 

This singularities can be modeled in a two dimensional 8-noded isoparametric 
element by 

* shifting the mid side nodes near the crack tip to the quarter point as shown in 
Fig. 4.1.b [16-19]. 
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* Collapsing an edge of the element to a single node at the crack tip [16-19] apart 
from, shifting the mid side nodes to the quarter point to form triangular crack tip 
elements as shown in Fig. 4.I.C. 

To model l/-\/r singularity the collapsed nodes are constrained to move as a 
single point and no relative movement between the nodes is allowed. To model 1 /r 
singularity three nodes which are constrained to move as a single point are released 
so that relative movement between them is possible. 

Using the above procedure the singularity is modeled in the interior as well as on 
the boundary of the element. 


4.2 Estimation of Fracture parameters 


In the study of structural problems involving cracks, the aim is to find the load at 
which the crack growth starts initiation leading to the ultimate failure of the structure. 
In EPFM J-integral and CTOD are the two widely used fracture parameters. 

The fracture parameter called the J-integral developed in [33] which is essentially 
the energy released per unit crack extension per unit thickness of the cracked body. 
This parameter is applicable to elastic or elasto-plastic behavior. 

The two dimensional form of the J-integral, with reference to crack problems, is 
given by 


J = i IWdy-T^^ds] 

with lU(e) = /q (7ij dcij where F is an arbitrary contour around the crack tip starting 
from one face and ending at the opposite face of the crack (fig. 4.2) 

ds is an element of arc along the integration contour, 

T is the traction vector. 

u is the displacement vector. 

The value of J is independent of the actual path chosen provided that the initial 
and end points of contour are on the opposite faces of the crack and that crack faces 
are stress free. 

It has also been shown [33] that J-integral, when defined along a contour around 


a.Parciit clt'jiu'iit 




c.Triangular quarter point element 


Figure 4.1: Cracktip singular elements 
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the crack tip represents the change in potential energy for unit crack area extension 
da. 


i.e J = 


9V_ 

da 


where V is the potential energy. Crack growth or fracture is said to occur for a 
particular load when the J value for the given loading conditions exceeds a critical 
value J\o which is a material constant. 


The method used for evaluation of Energy ReleaseRate in the present study is 
given below. 

Consider the load-displacement diagram for a crack of size ’a’ represented by the 
curve OA in Fig. 4.3. For a crack of size ’a-j-da’ the load displacement diagram is 
shown by the curve OE. Suppose the crack extension from a to a + da takes place 
for a particular load P. Then, the Energy Release Rate for that particular load P by 
the crack extension da is given by the area between the two curves divided by B da 
where B is the thickness of the body. 

Once the displacement field near the crack tip is known the crack tip opening 
displacement can be found out. 
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Chapter 5 


RESULTS AND DISCUSSION 


Some sample problems are solved using the program developed. They are discussed in 
this chapter. Three problems are solved for bodies having cracks. In these problems 
the effect of geometric nonlinearity on fracture parameters is reported. 


a) Cantilever Subjected to Moment at its Free End 


A geometric nonlinear analysis is carried out for a cantilever subjected to a moment 
at its free end. Fig. 5.1 gives the details of the geometry and the material properties 
used in the analysis. 

The moment is simulated using two equal and opposite forces at the free end of the 
cantilever as shown in Fig. 5.2. Since the loading is applied in steps the forces have to 
be follower forces. This is simulated in the program. The cantilever is modelled using 
20 8-noded isoparametric plane stress elements. The mesh is shown in Fig. 5.3. The 
same problem is analytically solved in [34]. Results obtained using our program are 
compared with the analytical solution. Good agreement with the results is evident 
from Table 5.1 and Table 5.2. Displacement is in cm and Moment is in N cm. 


Table 5.1: Comparision of Axial displacement 
Moment Analytical Present Work 


36.65 .6 .6003 

73.30 1.187 1.187 

109.95 2.06 2.06 

146.60 3.17 3.175 

183.25 4.8 4.45 
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Table 5.2: Comparision of Transverse displacement 


Moment 

Analytical 

Present Work 

36.65 

5.14 

5.14 

73.30 

10.15 

10.15 

109.95 

14.93 

14.93 

146.60 

20.00 

19.375 

183.25 

25.0 

23.45 


a) Cantilever Subjected to End Moment 




Dimensions 1 = 100 cm b = 20 cm h = 1 cm 
Properties E = 2.1 x 10'* Njcm^ i/ = 0.0 


Figure 5.1: Cantilever subjected to end moment 


Fig. 5.4 shows the load deflection response and Fig. 5.5 shows the deformed shape 
of the cantilever. 


b) Large deflection analysis of an Arch 

A large deflection analysis of an arch fixed at both ends is carried out using the 
program developed. The arch is subjected to a point load at the middle as shown in 
Fig. 5.6. Considering the symmetry of the problem only one half {left half) of the 
arch is modelled. The finite element mesh is composed of 12 8-noded isoparametric 
plane stress elements as shown in Fig. 5.7. Fig. 5.8 shows the final deformed shape 
of the arch. Fig. 5.9 shows the load deflection response of the arch. It can be seen 
that the arch stiffens with increasing load. 
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Figure 5.2: Moiuent on Cantilever 

Problems involving Cracked Bodies 


The Elastic Plastic material matrix whose formulation was explained in section 2.1.4 
was validated with the help of results available in literature. After such validation 
the effect of geometric nonlinearity on fracture paramters was investigated. 


c) Simply Supported Beam with a Central Crack 
under Point load 

A simply supported beam with a central crack subjected to central point load is 
investigated to study the effect of geometric nonlinearity on fracture parameters. 
Fig. 5.10 gives the dimensions and material properties used in the analysis. Po is the 
total applied load. It is assumed that the material is elastic and perfectly plastic. 

Fig. 5.11 shows the mesh used for finite element analysis. Considering the symme- 
try of the problem only one half of the beam is modelled. The structure is decornposed 
into 96 8-noded plane stress elements for analysis. To model crack tip singularity the 
analysis uses triangular quarter point 8-noded isoparametric elements as explained in 
section 4.1. Before adopting this mesh to study the effect of geometric nonlinearity 
on fracture parameters, a convergence study was carried out using three meshes. The 
first mesh had 64 elements and 231 nodes, the second mesh consisted of 96 elements 
and 335 nodes and the third mesh had 144 elements and 491 nodes. We observed a 



Figure 5.3: Finite element mesh of cantilever 
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Figure 5.6: Arch under Point load 
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Figure 5.8: Deformed shape of Arch 





c) Simply supported beam with a central 
crack under point load 

r = A/V. n=M5 



jL X 10^ N/cm^ 

(Ty — 200 Njcm? I 



Properties: E 




0.3 
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monotonic convergence and there was no mesh dependence on the results. Consid- 
ering computational efficiency and accuracy we adopted the second mesh i.e with 96 
elements and 335 nodes for further analysis. 

The effect of geometric nonlinearity on energy release rate is shown in Fig. 5.12 
and the effect of geometric nonlinearity on CTOD is shown in Fig. 5.13. The plastic 
zone near the crack is shown in Figs. 5.14a, 5.14b and 5.14c. The plastic zone at 
limit load is shown in Fig. 5.15 

d) Center cracked plate under uniaxial loading 


A problem of center cracked plate of finite size subjected to uniaxial loading is consid- 
ered here for analysis. The objective is to evaluate the effect of geometric nonlinearity 
on fracture parameters. 

The geometry and material properties of the plate are shown in Fig. 5.16. Po is 
the total applied load. Because of symmetry, a quarter of the plate shown in Fig. 5.16 
is taken for analysis. The finite element mesh is shown in Fig 5.17. The singularity 
at the crack tip is modelled by the use of triangular crack tip elements. It is assumed 
that the material is elastic and perfectly plastic. 

The effect of geometric nonlinearity on energy release rate and CTOD is shown 
in Figs. 5.18 and 5.19. The plastic zone near the crack is shown in Figs. 5.20a, 5.20b 
and 5.20c. 


e) Edge cracked plate under uniaxial loading 


A problem of edge cracked plate of finite size subjected to uniaxial loading considered 
here for analysis. The objective is to evaluate the effect of geometric nonlinearity on 
fracture parameters. 

The geometry and material properties of the plate are shown in Fig. 5.21. Po is 
the total applied load. Because of symmetry, half of the plate shown in Fig. 5.21 is 
taken for analysis. The finite element mesh is shown in Fig 5.22. The singularity at 
the crack tip is modelled by the use of triangular crack tip elements. It is assumed 
that the material is elastic and perfectly plastic. 

The effect of geometric nonlinearity on energy release rate and CTOD is shown 
in Figs. 5.23 and 5.24. The plastic zone near the crack is shown in Figs. 5.25a, 5.25b 
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Figure 5.11: Finite element mesh for the cracked beam 
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Figure 5.15: Spread of Plastic zone for cracked beam at limit load 







Figure 5.17: Finite element mesh for the center-cracked plate 



0 9 

Load N/cm 
Material : Elastic Per-Fectly Plastic 


Figure 5.18: Energy Release Rate vs Load for center-cracked plate 
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Figure 5.22: Finite element mesh for the Edge-cracked plate 
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and 5.25c. 


Cracked Bodies with Elastic Linearly Strain Hard- 
ening Material 

The problems of simply supported beam with a central crack, plate with center crack 
and plate with edge crack are solved in this section assuming that the material is 
elastic and linearly strain hardening in nature. The strain hardening co-efficient is 
taken to be H' = 2X10^. 

The effect of geometric nonlinearity on energy release rate and CTOD for simply 
supported beam with central crack is shown in Figs. 5.26 and 5.27. The plastic zone 
near the crack is shown in Figs. 5.28a, 5.28b and 5.28c. 

The effect of geometric nonlinearity on energy release rate and CTOD for plate 
with center crack is shown in Figs. 5.29 and 5.30. The plastic zone near the crack is 
shown in Figs. 5.31a, 5.31b and 5.31c. 

The effect of geometric nonlinearity on energy release rate and CTOD for plate 
with edge crack is shown in Figs. 5.32 and 5.33. The plastic zone near the crack is 
shown in Figs. 5.34a, 5.34b and 5.34c. 

It is observed that for same load the body with elastic perfectly plastic material 
has higher energy release rate and larger CTOD than the body with elastic linearly 
strain hardening material. 



Lo ad N/cm 


Material : Elastic Perfectly Plastic 


Figure 5.23: Energy Release Rate vs Load for Edge-cracked plate 







a) P =z 16.2 


6) P = 32.4 



c) P = 40.5 


Figure 5.25: Spread of Plastic zone 


for Edge-cracked plate 







6] 




a) P = 57.75 


b) P = 80.85 



Figure 5.28: Spread of Plastic zone for cracked beam 





Load Njcm 


Material : Elastic Linearly Strain Hardening 


Figure 5.30: CTOD vs Load for center-cracked plate 
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Lo ad N/cm 


Material : Elastic Linearly Strain Hardening 


Figure 5.33: CTOD vs Load for edge-cracked plate 






c) p = 48.6 


Figure 5.34: Spread of Plastic zone for edge-cracked plate 


Chapter 6 


CONCLUSIONS AND SCOPE 
FOR FUTURE WORK 


The following conclusions can be drawn from the present study 

• Geometric and Material nonlinearities have been successfully modelled. 

• The nonlinearity due to large deformations in the structure influences the load 
displacement characteristic to a great extent. 

• Geometric nonlinearity has influence on the fracture parameters. Therefore to 
have a more accurate design Geometric nonlinearity should be considered in 
evaluating fracture parameters. 

Scope For Future Work 

• Modified Newton raphson technique takes significant time as we go for mesh 
refinement in our analysis. Inclusion of some acceleration scheme in this method 
would make the program computationally more efficient. 

• Present work can be extended to include the dynamic analysis of cracked struc- 
tures thus making it more general. 
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Appendix 

Matrices used in two dimensional analysis in T.L formulation 

Incremental Strains 


0^11 = oWl,l + oUi,i + 0^2,1 0*^2, 1 + 2 ((0^1,1)^ + (o«2,l)^) 

0^22 = 0^2,2 + 0^1,2 0^1,2 + o^2,2 0t^2,2 + 2 ((o^l.s)^ + (o«2,2)^) 

0^12 = 2 + 0«2,l) + 2 (o“l,l 0 Mi,2 +0 ^2,1 0^2,2 +0 ^1.2 0«1,1 +0 ^2,2 0^2,1) 

+ 2 ( 0 ^ 1,1 0 ^ 1,2 + 0 ^ 2,1 01 ^ 2 , 2 ) 

where o«i,y = = ftj; 

Linear Strain-Displacement Transformation Matrix 

using oe = IBl u 

where oe^ = [oen 0622 2oei2] 

= [nl U 2 “2 • • • '^2] 

qBl = qBlo + o^Ll 

o^jv,i 0 
0 0^7V,2 

ohjv,2 O^iV,! 


^11 0^1,1 hi 0^1,1 • ■ • hi O^AT,! 

I 12 0 ^ 1,2 ^22 0 ^ 1,2 ••• ^22 0 ^JV ,2 

Al A2 ... A3 



IBlo = 


0^1,1 0 0^2,1 0 0^3,1 0 

0 o/ii ,2 0 0 ^ 2,2 0 0 ^ 3,2 

[ 0^1,2 ohi,i 0^2,2 0^2,1 0^3,2 0^3,1 


where o/ifcj = - ‘w • ; 


^^1 — Y^k=:i N = Number of nodes; 
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where 


A 1 — {111 0^1, 2 + ^12 ohi,i) A 2 = {I21 0^1, 2 + ^22 0^1,1) -A 3 = (/21 o^jv, 2 + ^22 o^iv,i) 
where 

N N 

hl = ^ 2 ohk,l *«!,■ l 22 = Y]ohk ,2 

k=l k=l 

N N 

h\ = '^ohk,i *^ 2 ] h2 = '^ohk,2 

k=l k=l 

N = Number of nodes; 

Nonlinear Strain-Displacement Transformation matrix 


qBnL = 


0^1,1 0 0^2,1 0 0^3,1 0 • • • O^JV.l 

0^1,2 0 0^2,2 0 0^3,2 0 ... ohN,2 

0 0 ^ 1,1 0 0 ^ 2,1 0 0 ^ 3,1 • • • 0 

0 0^1,2 0 0^2,2 0 0^3,2 • • • 0 


0 

0 

O^N,! 

ohN,2 


Second Piola-Kirchoff Stress Matrix and Vector 


'■‘S'n iS: 


t c _ 
0-0 — 


0 0 

0 0 

0 0 o5ii qSi2 

0 0 q5'21 qS22 


0*^11 0*^12 
IS 21 IS 22 


t r^C ^ C t C ^ 

qO — [0^11 0^22 0*^12 J 

Determination of the ratio of mass densities 

The mass density ratio °/3/V can be evaluated since the mass of the particles 
considered is conserved: 


[ V (i*xi ctx2 d*X3 = / cPxi <fx2 dPxz 
J*v J°v 


But 


d^Xi d*X2 d*X3 = {det qX) dPxi dPx2 d^xa 

T 

where qX is the deformation gradient given by (oV ^X'^) 


73 



where 

Hence, 


T _ r d 


oV^ = [ 


a 


d^xi d^X2 


d 

d^xz 


] and = [ *Xi ^X2 ^xz 


V = V {dttlX) 
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